Edge currents in superconductors with a broken time— reversal symmetry 
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We analyze edge currents and edge bands at the surface of a time-reversal symmetry breaking 
d x 2_ y 2 + id xy superconductor. We show that the currents have large Friedel oscillations with two 
interfering frequencies: \/2kF from sub-gap states, and 2Uf from the continuum. The results are 
based independently on a self-consistent slave-boson mean field theory for the t — J model on a 
triangular lattice, and on a T-matrix scattering theory calculation. The shape of the edge-state 
band, as well as the particular frequency V2kp of the Friedel oscillations are attributes unique for 
the d x 2_ y 2 + id xy case, and may be used as a fingerprint for its identification. Extensions to different 
time-reversal symmetry breaking superconductors can be achieved within the same approach. 



Superconductors which break time-reversal symmetry 
have attracted a great deal of attention recently. The 
most prominent example is S^RuOzi, which has a p- type 
complex order parameter Q. There have been theoreti- 
cal suggestions that Na^CoCVyH^O belongs to this class, 
but with complex d-wave symmetry 0, 0, ■ It is known 
that interesting edge effects are induced by the surface 
in such superconductors. First, the surface induces the 
appearance of an edge band inside the superconducting 
gap. Second, a chiral edge current parallel to the surface 
appears. As shown in several independent semi-classical 
calculations HHI1||, the shape of the dispersion e& 
of the edge states depends on the symmetry of the su- 
perconductor and may be used as a fingerprint for the 
latter. In this letter, we quantitatively address this prob- 
lem and point out a novel interference in the edge cur- 
rent which distinguishes between d x 2_ y 2 + id xy (= d + id') 
and complex p time-reversal symmetry breaking super- 
conductors. We perform a self-consistent solution of the 
t — J model on a triangular lattice which has a d + id' 
superconducting phase as the solution of its slave-boson 
mean- field theory 0, 0- EI ■ In order to explain the novel 
results uncovered by the numerical solution, we analyti- 
cally solve a quantum continuum model by treating the 
quasiparticle scattering on the surface using an extension 
of the T-matrix formulation of 

The results are summarized as follows (see Fig.^l: The 
surface induces edge currents and an edge band in agree- 
ment with the semi-classical prediction of two parts of a 
parabola et(k x ) = —Asgii(k x )(2k x /k F ~ 1). As a sur- 
prising novel result, the size and direction of the current 
changes rapidly with distance from the surface. This is 
interpreted as Friedel oscillations of two frequencies, 2kp 
and \/2kp. The 2k p oscillations are the usual Friedel 
oscillations of the continuum states. The \/2kp oscilla- 
tions are due to the zero energy mode of the sub-gap 
band, ej, = 0, and correspond to a quasiparticle scatter- 
ing on the surface at an incident angle of 45° [g] (Fig. 
(a) inset). Most notably, the amplitudes of both oscilla- 
tions are identical, and comparable to the nonoscillating 



part of the current; from the self-consistent numerical re- 
sults we see that they are indeed large enough to reverse 
the current direction in some regions close to the surface. 
The overall magnitude of the edge current, therefore, is 
considerably smaller than could be expected naively. Our 
calculations lead to an integral edge current I of the order 
of 60nA, which also corresponds roughly to its maximal 
amplitude. Such a current may be captured by measur- 
ing the induced magnetic field B. For our model cal- 
culations, we estimate its maximum at the surface as 
B = /.iqI/2tt£ ~ 0.1G, where £ is the coherence length of 
the superconductor. 

The t — J model on a two-dimensional triangular lat- 
tice is an example of a microscopic model with a d + id' 
superconducting phase as the result of the slave-boson 
mean-field theory [3- 0- Q • We solve the mean-field the- 
ory self-consistcntly in real space following the approach 
by [Hi m the presence of two surfaces on the top and 
bottom of the system. The model Hamiltonian is the 
t — J model plus a long-range Coulomb repulsion, 
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where (ij) runs over the nearest neighbors lattice sites 
of an equilateral triangular lattice, Ci„ are the electron 
operators, and Si the spin 1/2 operator. The model 
is completed by the constraint of no double occupancy 
per site, n, = J2a c i<j c i<? < 1- J > is the antiferro- 
magnetic exchange interaction, and t the hopping inte- 
gral. Because of the missing particle-hole symmetry on 
the triangular lattice, the sign of t is important 0, 0- 
We choose t/J = —3, corresponding to electron dop- 
ing of the system. V t = V^j^Mj ~ ~ *j\> is 
the Coulomb potential felt by the particle at the site i, 
where n is the average density, j runs over all lattice sites, 
and V c s» 5 J ^(J- The long-range Coulomb interaction is 
necessary to overrule the inherent tendency of the mean- 
field theory to phase separate. In a uniform system, Vi 
vanishes. On the lightly doped triangular lattice, how- 
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ever, the Coulomb interaction alone is not strong enough 
against the phase separating instabilities. The key in- 
sight for a stabilization is the correction of the mean-held 
theory by a Jastrow-like modification of the hopping pa- 
rameter t — > texp(—w\y/(l — rii){l — rij) — (1 — n)\) for 
nearest neighbor sites with w ~ 5. This factor suppresses 
the gain in kinetic energy by clustering carriers on neigh- 
boring sites. For small fluctuations in the carrier concen- 
trations, it is close to 1. Its precise form otherwise has 
not much influence. 

In the slave-boson formulation, the electron operator 
is decomposed as c ia = fi a b\, where fi„ is a spin carrying 
fermion, and bi a charge carrying boson operator. The 
constraint rij < 1 becomes the identity ^ CT f} a fia+b\bi = 
1, which must be fulfilled at each site i, and which can 
be included in the action with Lagrange multipliers A,;. 

In the superconducting phase, the bosons bi are con- 
densed and directly related to the local carrier concen- 
trations, Xi = (1 — rii) — b\. The (Jastrow-corrected) 
mean-field theory can be derived by a variational Ansatz 
(llj . and leads to a Hamiltonian for the fermion fields 
fi a , expressed in terms of the X4, the hopping param- 
eter Xij — (fhfje)' an< l the pairing parameter = 
(filfjl ~ filfjf)- The theory is solved self-consistently 
by iteration: For random initial order parameters, we 
solve the Bogoliubov-de Gennes equations by diagonal- 
ization of the Hamiltonian. The resulting spectrum and 
wave functions are used to reconstruct the order parame- 
ters and the local densities. The Lagrange multipliers Xi 
are calculated explicitly by minimizing the action with 
respect to the condensed bosons bi. This is repeated un- 
til convergence and requires for the triangular lattice at 
low dopings (< 15%) about 500 to 1000 iterations. 

We study lattices consisting in a strip with periodic 
boundary conditions along the x-direction and, to model 
the two surfaces, open boundary conditions along the y- 
direction. To reduce the interference between the edges, 
we choose systems in which the surfaces are far apart 
(typically ^100 lattice sites). Finite size constraints are 
important and distort the perfect d + id'— symmetry in 
that the phase difference of Ajj between lattice directions 
deviates from 2ir/3, or that the ratio ReA^ /ImA^ in the 
bulk deviates from unity. To minimize these effects, we 
choose a sufficiently large extension along the x-direction 
of N x = 24. The computational expense then limits the 
y-extension to about 100 sites. 

In Fig. ^ ( a ) (circles) we show the resulting edge cur- 
rent. It has oscillations with a large amplitude which 
lead at some layers to the reversal of the current direc- 
tion. The spatial period of the oscillations is roughly 
twice the lattice spacing, but beatings indicate a fre- 
quency mixture. The multiplication of the displayed cur- 
rent by te/h (e = electron charge) provides the current in 
Amperes. With t = —3 J and J = 20mcV (as estimated 
for the Na^CoO^ superconductor, which was considered 
as a candidate for the d+id' symmetry 0,0,0|)i the over- 



all current, integrated along the direction perpendicular 
to the surface, is about 60nA (setting an order of magni- 
tude), and is concentrated within a few times the decay 
length of £ cur « 5a (with a w 3A the lattice constant for 
the equilateral triangular lattice). 

Part (b) of the figure shows the spectrum of the system 
as a function of k x after a Fourier transformation of the 
real-space result. We have suppressed the edge states at 
one of the two surfaces by inspection of the support of the 
wave functions. An asymmetric sub-gap band is clearly 
visible. Its shape deviates slightly from the parabola pre- 
dicted from the semi-classical arguments as an effect of 
the discrete and finite lattice and the variation of the 
magnitude of A close to the edge. In particular, a small 
gap persists as the result of the limited system size and 
the interference between the edges. The gap closes with 
increasing distance of the surfaces. 

A similar computation can be performed for a lattice 
with a single site defect 01 Instead of a band, a pair 
of bound states appears at the defect site. For similar 
interaction parameters, the edge current is concentrated 
on the immediate neighbors of the defect site only, and 
has a value of about lOOnA. This leads to a magnetic 
field of the order of 1G on the defect site. 

The large amplitude of the oscillations seems to be sur- 
prising, but can be derived within a continuum model for 
the scattering on the surface. We show that the oscilla- 
tions are Friedel oscillations induced by the scattering on 
the surface. The calculation must go beyond the semi- 
classical arguments, and we follow the T-matrix calcu- 
lation by |9( , which we extend for the calculation of the 
current density j(y). Even though we focus on the d+id' 
symmetry, the approach is general and can be extended 
to any given symmetry of the superconductor. 

We consider a semi-infinite two-dimensional d+id' su- 
perconductor in the (x, y)-plane with its surface at y = 0, 
and the superconductor at y > 0. The current can be ob- 
tained from the causal (time-ordered) Green's function, 
G c (x, x';t), by the expression 

J'W = ^E( V * - V X ')<^(x,x';t = 0-)| x , =x , (2) 

a 

with a the spin projections, and m the mass of the carri- 
ers. In contrast to causal instead of retarded Green's 
functions must be used. 

For the d + id' superconductor, we write the gap func- 
tion in the simplified, position independent form Ak = 
Ae 2 * e , where 9 is the angle of the momentum k to a given 
axis. Since we intend here to provide an explanation to 
the observed oscillations while keeping the calculation as 
simple as possible, we assume |A| to be constant even in 
the vicinity of the surface. The variability of A is taken 
into account, however, in the self-consistent numerical 
solution of the Bogoliubov-de Gennes equations. Since 
the amplitude of the gap is independent of the direction, 
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the orientation of the surface is of no importance. In 
this approximation, the Fermi surface is assumed to be 
cylindrical. The surface is modeled by a line of scatter- 
ing centers with elastic potentials V, which eventually we 
let tend to infinity. The preserved translation symmetry 
along the x-axis allows us to keep the momenta parallel 
to the surface, k x , so that the scattering equation reads 

Gl (y, y>; uj) = (y-y>; uj)+gl (y; uj)T kx (w)g^ (-</; uj), 

_ (3) 

where the bold symbols denote functions in the particle- 
hole Nambu space, and uj is the frequency. The free 
causal Green's function, g£ (y), is 



nN 



-e 4 ?( 



pVA 2 - uj 2 - iO 
(cjcos(py)xo + A ^2 cos(py + v26)t v }, (4) 



where N is the density of states at the Fermi energy, kp 
the Fermi vector, p = yfcj — k 2 , and £ the coherence 
length. Here k x = kp cos(6),p — kpsm(9), 8 S [0, ir]. 
The t± — t x ± iT y are the Pauli matrices in the Nambu 
space; To is the unit matrix. 

The T-matrix is defined by T^(u)) = Vt z - g^(o;). 
To model the surface, we let V — > oo. With £1 = 



VA 2 



iO and Eq. |0J this leads to 



pil/irNo 



A 2 cos 2 (26l)(l -»0) V-Acos(26>) 



-Acos(20) 
(5) 

In this expression we still have scattering between phys- 
ical y > states and unphysical y < states. It con- 
tains two sub-gap bands with the energies ' (k x ) — 
±Acos(26 l ) for the states defined for y > and y < 
respectively. The T-matrix must be split into two parts 
acting on those states separately. For the time-reversal 
symmetry breaking superconductor, we expect an asym- 
metry of these bands with respect to the sign of k x , lead- 
ing to the only possible choice 

c-b\kx) = -sgn(fc*)Acos(20) = -sgn(fc x )A(2^ - 1), 

(6) 

This is precisely the expression ob- 
tained from the semi-classical models 0, § , and we can 
focus on the sub-gap band e^\k x ) = £b(k x ) only, rel- 
evant for y > 0. The splitting is achieved by writ- 
ing l/[(w — )(u> — e& )] as the sum and difference 

( 1 2) 

of l/(w — ), such that the remaining factors in the 
T-matrix (JSJ) compensate any unphysical singularity ap- 
pearing from this decomposition. The T-matrix be- 
comes, disregarding the y < part, 



and £j = — 



pn/2nN ( 1 sgn(fc,; 
uj - e b (k x )(l - iO) \sgn(k x ) 1 



(7) 



In (k x , y, w)-space, the current density (0) then reads (us- 
ing Nq = m/2ir) 



2ik x 



sL (v; w)T*. (w)g£ (-y;u))]u 



j(k x ,y;u) ■■ 

~ 2pft w - e 6 (l - iO) 
- Asin(2p2/) |sin(20)| (w - e 6 ) + A 2 sin 2 ( ra ) sin 2 (26»)|. 



|cos 2 (py)(w - e 6 ) 5 



(8) 

The w-integration leads to two contributions, the sub- 
gap poles, which exist for uj = t\, < only, and an integral 
running along the continuum state branch cut at uj < 
—A. For the sub-gap states, we obtain (see also @): 



f s (y) = -7- 



A ^ dfc > e -f|cosW|„,,i 



4tt 



— kf 



sin (py) 
x |sin(20)|e(-e 6 ), (9) 



with 9 the unit step function. Close to the surface, y <C 
£, we can replace the exponential by unity. At y w 0, 
the current grows as y 2 . For yfcf ^> 1 but still y <C £, 
the integrand oscillates rapidly, and, to leading order in 
1/kpy, j sg can be estimated as 



f g (y) 



Ak f 



4tt 



1 



i(V2k F y) 



3%/2 



2k F y 



(10) 



i.e. the spatial oscillations are determined by the zero 
energy mode of quasiparticle scattering on the surface, 
for which p = kp/\2 |6( (see Fig.Q](a), inset). Constant 
contributions at the artificial cutoffs at k x = ±kp have 
been neglected. For y exceeding the coherence length £, 
this current is suppressed exponentially. This allows us 
to identify £ cur = £. 

In a time-reversal symmetry breaking superconductor 
the continuum states carry a current, too. Their con- 
tribution is given from the w-integral running along the 
branch cut at (— oo,— A). We restrict the integration 
over k x to positive values, and keep only terms which are 
even in k x . This yields (uj — > — uj) 



j (y) = — 



dk x 

Kx^b 

p 



dw 



c°s(ff A /(IF^l) 



A 



Vw 2 - A 2 



(uj 2 - A 2 ) cos 2 (py) + A 2 sin 2 (26l) cos(2py) 
(uj 2 - A 2 ) + A 2 sin 2 (26>) ' 



(11) 



The integral weight is dominated by the singularity 
\I\Juj 2 — A 2 and is concentrated at the physically rel- 
evant values of uj rj A. In this region, the last factor in 
Eq. l(TT|l varies slowly about cos(2py), and we replace it 
by this quantity. For the w-integration provides a 

Bessel function which accounts for the required exponen- 
tial decay. At y = 0, however, the w-integration reduces 
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to the integration of l/y/u> 2 — A 2 , which is divergent. 
This is an artifact of the approximation of a constant No. 
The integral must be cut off at the band-edge where the 
density of states vanishes and, therefore, is given solely 
by its value at lu — A, J dui/ \Ju 2 — A 2 |a = tt/2. For 
y close to the surface, we replace the w-integral by this 
value. The remaining k x integration is elementary. We 
obtain for y <C £ , to leading order in \jkpy, 

Ak F sm(2k F y) 

3 (y) = -j 777— — • 12 

An 2k F y 

The full current density, j(y), is the sum of Eqs. (|10|) and 
(|12[1 . and shows Friedel-like spatial oscillations with the 
two frequencies \2kp and 2k F . 

Even though details vary in real physical systems and, 
as shown above, in fully self-consistent solutions, these 
expressions make the following precise statements: i) The 
edge state dispersion is identical to the semi-classical re- 
sults 0,0. ii) There are two interfering Friedel oscilla- 
tions of the same amplitude with the frequencies \2kp 
and 2k F , where the former is determined by the condi- 
tion Cb(k x ) — 0. Hi) Most notably, the amplitude of the 
oscillations is comparable to the non-oscillating term in 
the current, i.e. it can reverse the direction of the cur- 
rent in some regions close to the surface. The exact ratio 
of the amplitudes cannot be captured with the present 
(not self-consistent) calculation. Yet the self-consistent 
numerical calculation above shows that this is indeed the 
case (Fig. iv) The current amplitude is proportional 
to A and kp. 

Note that Friedel oscillations with the same two fre- 
quencies are found for the carrier density, n(y), since 
it connects to the Green's function through n(y) = 
-2i[G%,t/;t = 0-)]n 

We also remark that in the spin-zero sub-gap band of 
a p x + ip y superconductor with Afc = Ae* e |9| the same 
calculation leads to £b ~ k x , which vanishes at p = kp. 
This provides a frequency 2k p for the Friedel oscillations 
from the sub-gap contribution and is indistinguishable 
from the continuum state oscillations. 

From the spectrum in Fig.Q](b), we obtain k x w 2.2/a 
at = A, and k x w 1.6/a at eb — 0, with a ratio of pre- 
cisely \J~2~. With these two values, and A taken from the 
bulk of the lattice, we can superpose the numerical result 
with the theoretical expressions and p2|l . We com- 
plete the theoretical expression by an exponential decay 



factor e~ y ^ cur . While the the amplitude and frequencies 
are fixed by A and kp, we fit the data for £ cur , the ef- 
fective distance between the surface and the first lattice 
layer, and the value of the nonoscillating part (replacing 
the l/3\/2). The resulting curve shows a nice agreement 
between the theoretical and numerical currents (Fig. Q] 
(a)). We deduce a £ cur « 5a, and a nonoscillating part, 
which is an order of magnitude below the prediction of 
1/3 v2- The integration of both the theoretical and nu- 
merical currents leads to an integral current of 60nA. 

To conclude, we have shown with two independent 
calculations that two types of Friedel oscillations exist 
in the time-reversal symmetry breaking superconductor, 
and that they play an essential role in the surface effects. 
The expected integral edge current and the induced mag- 
netic field is reduced by about an order of magnitude with 
respect to the estimates for a nonoscillating prediction, 
which has important consequences for experiments. 
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FIG. 1: (a) Current (in dimensionless units) from the 
self-consistent calculation (circles) and the prediction of 
Eqs. IIUI and 1121 multiplied by an exponential decay factor 
exp(— y/£cur) (full line). The inset shows the quasiclassical 
Andreev-scattering on the surface leading to the \/2&f oscil- 
lations (full lines: particle scattering, dashed lines: hole scat- 
tering), (b) Spectrum -E n> fc x as a function of the momentum 
k x \ n is the quantization along the y-direction. Continuum 
states form the hat-like structures, the sub-gap states have 
roughly parabolic forms. The small gap at zero energy is a 
finite size effect due to the coupling between the two edges. 



